In this script we are going to compute marker genes for the samples at varying resolutions.
library(Seurat)
library(ggpubr)
library(cowplot)
library(dplyr)
library(ggplot2)
library(stringr)
library(readr)
Loading necessary paths and parameters
source(here::here("misc/paths.R"))
source(here::here("utils/bin.R"))
set.seed(123)
source(here::here("misc/paths.R"))
source(here::here("utils/bin.R"))
"{anot}/{plt_dir}" %>%
glue::glue() %>%
here::here() %>%
dir.create(path = .,
showWarnings = FALSE,
recursive = TRUE)
"{anot}/{robj_dir}" %>%
glue::glue() %>%
here::here() %>%
dir.create(path = .,
showWarnings = FALSE,
recursive = TRUE)
The data used in this Rmarkdown document comes from 03-clustering_integration.Rmd where the data was integrated.
merged_se <- "{clust}/{robj_dir}/integrated_spatial.rds" %>%
glue::glue() %>%
here::here() %>%
readRDS(file = .)
First we are going to look at the cluster resolution we are going to check
Seurat::DimPlot(merged_se, group.by = "Spatial_snn_res.0.3")
res_vec <- c("Spatial_snn_res.0.01", "Spatial_snn_res.0.05",
"Spatial_snn_res.0.1", "Spatial_snn_res.0.3",
"Spatial_snn_res.0.8", "Spatial_snn_res.1",
"Spatial_snn_res.1.25", "Spatial_snn_res.1.5")
# Iterate over resolutions
lapply(res_vec, function(i) {
# Iterate over slices
lapply(id_sp_df$gem_id, function(id) {
Seurat::SpatialPlot(
merged_se,
group.by = i,
images = id,
image.alpha = 0,
crop = FALSE,
pt.size.factor = 1.25)
}) %>%
patchwork::wrap_plots(., nrow = 2, guides = "collect")
})
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
We are going to use the function FindAllMarkers to get the marker genes for all the clusters at the varying resolutions:
# Assign Identities
Seurat::Idents(merged_se) <- merged_se@meta.data[, "Spatial_snn_res.0.3"]
# Compute markers
markers03 <- Seurat::FindAllMarkers(
object = merged_se,
assay = "Spatial",
slot = "data",
only.pos = TRUE)
DT::datatable(
data = markers03,
filter = "top")
out_file03 <- "{anot}/{robj_dir}/markers_Spatial_snn_res.0.3.xlsx" %>%
glue::glue() %>%
here::here()
# Remove file to write it from scratch
file.remove(out_file03)
lapply(unique(as.character(markers03$cluster)), function(clust) {
markers03 %>%
dplyr::filter(cluster == clust) %>%
dplyr::arrange(dplyr::desc(avg_log2FC)) %>%
xlsx::write.xlsx(
.,
file = out_file03,
row.names = FALSE,
sheetName = glue::glue("Cluster-{clust}"),
append = file.exists(out_file03))
})
We are going to use the function FindAllMarkers to get the marker genes for all the clusters at the varying resolutions:
# Assign Identities
Seurat::Idents(merged_se) <- merged_se@meta.data[, "Spatial_snn_res.1.5"]
# Compute markers
markers_15 <- Seurat::FindAllMarkers(
object = merged_se,
assay = "Spatial",
slot = "data",
only.pos = TRUE)
DT::datatable(
data = markers_15,
filter = "top")
out_file_15 <- "{anot}/{robj_dir}/markers_Spatial_snn_res.1.5.xlsx" %>%
glue::glue() %>%
here::here()
# Remove file to write it from scratch
file.remove(out_file_15)
lapply(unique(as.character(markers_15$cluster)), function(clust) {
markers_15 %>%
dplyr::filter(cluster == clust) %>%
dplyr::arrange(dplyr::desc(avg_log2FC)) %>%
xlsx::write.xlsx(
.,
file = out_file_15,
row.names = FALSE,
sheetName = glue::glue("Cluster-{clust}"),
append = file.exists(out_file_15))
})
merged_se[["annotation-general"]] <- dplyr::case_when(
merged_se@meta.data$Spatial_snn_res.0.3 == 0 ~ "T cell zone",
merged_se@meta.data$Spatial_snn_res.0.3 == 1 ~ "T:B Border",
merged_se@meta.data$Spatial_snn_res.0.3 == 2 ~ "Epithelial-1",
merged_se@meta.data$Spatial_snn_res.0.3 == 3 ~ "Inter-Follicular Zone",
merged_se@meta.data$Spatial_snn_res.0.3 == 4 ~ "Proliferating Follicle",
merged_se@meta.data$Spatial_snn_res.0.3 == 5 ~ "Germinal Center",
merged_se@meta.data$Spatial_snn_res.0.3 == 6 ~ "Subepithelial",
merged_se@meta.data$Spatial_snn_res.0.3 == 7 ~ "Epithelial-2",
merged_se@meta.data$Spatial_snn_res.0.3 == 8 ~ "Muscle"
)
Furthermore, we annotate these slides following a pathologists histological annotation. As of now we only have annotation of specific regions from one sample:
# load annotations
ann_dir <- "{anot}/{robj_dir}/mantle-zone-and-vascular-spaces" %>%
glue::glue() %>%
here::here()
hist_df <- lapply(list.files(ann_dir), function(n) {
id <- substr(n, 1, 13) # Extract first 13 characters
# Go into the directory for each slide
z_ls <- list.files(glue::glue("{ann_dir}/{n}"), full.names = TRUE)
tmp_df <- lapply(z_ls, function(z) {
# Iterate over region directories
fn_ls <- list.files(z, full.names = TRUE)
z_bc <- lapply(fn_ls, function(fn)
readr::read_csv(file = fn) %>% dplyr::pull("barcode")) %>%
unlist()
# Add region to data frame
dirs <- stringr::str_split(string = z, pattern = "/", simplify = TRUE)
zone <- dirs[, ncol(dirs)]
data.frame("barcodes" = z_bc, "histological_annotation" = zone)
}) %>% dplyr::bind_rows()
# Add gem_id to dataframe
tmp_df$gem_id <- id
tmp_df
}) %>%
dplyr::bind_rows()
# Make sure there are no duplicate rows
hist_df <- dplyr::distinct(hist_df)
Add histological annotation to seurat object
merged_se@meta.data <- merged_se@meta.data %>%
tibble::rownames_to_column("barcodes") %>%
dplyr::left_join(hist_df, by = c("barcodes", "gem_id")) %>%
dplyr::mutate(
histological_annotation = dplyr::if_else(is.na(histological_annotation),
"unannotated", histological_annotation)) %>%
tibble::column_to_rownames("barcodes")
table(merged_se@meta.data$histological_annotation)
##
## Mantle zones unannotated Vascular spaces
## 165 15971 88
Visualize the annotation
library(colorBlindness)
SpatialPlot(
object = merged_se,
group.by = "histological_annotation",
images = c("esvq52_nluss5", "exvyh1_66caqq"),
crop = FALSE, pt.size.factor = 1.25,
cols = c("#009292", "lightgrey", "#ff6db6"))
"{anot}/{robj_dir}/integrated_spatial_annot.rds" %>%
glue::glue() %>%
here::here() %>%
saveRDS(object = merged_se, file = .)
sessionInfo()
## R version 4.0.1 (2020-06-06)
## Platform: x86_64-conda_cos6-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux Server release 6.7 (Santiago)
##
## Matrix products: default
## BLAS/LAPACK: /scratch/groups/hheyn/software/anaconda3/envs/spatial_r/lib/libopenblasp-r0.3.12.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=es_ES.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=es_ES.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=es_ES.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] colorBlindness_0.1.9 readr_2.1.1 stringr_1.4.0 dplyr_1.0.7 cowplot_1.1.1 ggpubr_0.4.0 ggplot2_3.3.5 SeuratObject_4.0.4 Seurat_4.0.5
##
## loaded via a namespace (and not attached):
## [1] Rtsne_0.15 colorspace_2.0-2 ggsignif_0.6.3 deldir_1.0-6 ellipsis_0.3.2 ggridges_0.5.3 rprojroot_2.0.2 spatstat.data_2.1-0 farver_2.1.0 leiden_0.3.9 listenv_0.8.0 bit64_4.0.5 ggrepel_0.9.1 fansi_0.5.0 codetools_0.2-18 splines_4.0.1 knitr_1.36 polyclip_1.10-0 jsonlite_1.7.2 broom_0.7.10 ica_1.0-2 cluster_2.1.2 png_0.1-7 uwot_0.1.11 shiny_1.7.1 sctransform_0.3.2 spatstat.sparse_2.0-0 compiler_4.0.1 httr_1.4.2 backports_1.4.1 assertthat_0.2.1 Matrix_1.4-0 fastmap_1.1.0 lazyeval_0.2.2 cli_3.1.0 later_1.3.0 htmltools_0.5.2 tools_4.0.1 igraph_1.2.10 gtable_0.3.0 glue_1.5.1 RANN_2.6.1 reshape2_1.4.4 Rcpp_1.0.7 carData_3.0-4 scattermore_0.7 jquerylib_0.1.4 vctrs_0.3.8 nlme_3.1-153 lmtest_0.9-39 xfun_0.29 globals_0.14.0 mime_0.12 miniUI_0.1.1.1
## [55] lifecycle_1.0.1 irlba_2.3.5 rstatix_0.7.0 goftest_1.2-3 future_1.23.0 MASS_7.3-54 zoo_1.8-9 scales_1.1.1 vroom_1.5.7 spatstat.core_2.3-2 hms_1.1.1 promises_1.2.0.1 spatstat.utils_2.3-0 parallel_4.0.1 RColorBrewer_1.1-2 yaml_2.2.1 reticulate_1.22 pbapply_1.5-0 gridExtra_2.3 sass_0.4.0 rpart_4.1-15 stringi_1.7.6 highr_0.9 rlang_0.4.12 pkgconfig_2.0.3 matrixStats_0.61.0 evaluate_0.14 lattice_0.20-45 ROCR_1.0-11 purrr_0.3.4 tensor_1.5 labeling_0.4.2 patchwork_1.1.1 htmlwidgets_1.5.4 bit_4.0.4 tidyselect_1.1.1 here_1.0.1 parallelly_1.29.0 RcppAnnoy_0.0.19 plyr_1.8.6 magrittr_2.0.1 R6_2.5.1 generics_0.1.1 DBI_1.1.1 pillar_1.6.4 withr_2.4.3 mgcv_1.8-38 fitdistrplus_1.1-6 survival_3.2-13 abind_1.4-5 tibble_3.1.6 future.apply_1.8.1 car_3.0-12 crayon_1.4.2
## [109] KernSmooth_2.23-20 utf8_1.2.2 spatstat.geom_2.3-1 plotly_4.10.0 tzdb_0.2.0 rmarkdown_2.11 grid_4.0.1 data.table_1.14.2 digest_0.6.29 xtable_1.8-4 tidyr_1.1.4 httpuv_1.6.4 gridGraphics_0.5-1 munsell_0.5.0 viridisLite_0.4.0 bslib_0.3.1